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Abstract 

We characterize the maximum controlled invariant (MCI) set for discrete- as 
well as continuous-time nonlinear dynamical systems as the solution of an infinite- 
dimensional linear programming problem. For systems with polynomial dynamics 
and compact semialgebraic state and control constraints, we describe a hierarchy of 
finite-dimensional linear matrix inequality (LMI) relaxations whose optimal values 
converge to the volume of the MCI set; dual to these LMI relaxations are sum-of- 
squares (SOS) problems providing a converging sequence of outer approximations 
to the MCI set. The approach is simple and readily applicable in the sense that 
the approximations are the outcome of a single semidefinite program with no addi- 
tional input apart from the problem description. A number of numerical examples 
illustrate the approach. 

1 Introduction 

Given a controlled dynamical system described by a differential (continuous-time) or 
difference (discrete-time) equation, its maximum controlled invariant (MCI) set is the 
set of all initial states that can be kept within a given constraint set ad infinitum using 
admissible control inputs. This set goes by many other names in the literature, e.g., 
viability kernel in viability theory [5], or (A, £?)-invariant set in the linear case [14] . 

Set invariance is an ubiquitous and essential concept in dynamical systems theory, as far 
as both analysis and control synthesis is concerned. In particular, by its very definition, 
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the MCI set determines fundamental limitations of a given control system with respect 
to constraint satisfaction. In addition, there is a very tight link between invariant sets 
and (control) Lyapunov functions. Indeed, sub-level sets of a Lyapunov function give 
rise to invariant sets. Conversely, at least in the linear case, any controlled invariant set 
gives rise to a control Lyapunov function, and therefore these sets can be readily used to 
design stabilizing control laws; see, e.g., [9] for a general treatment and, e.g., [HI EZ] for 
applications in model predictive control design. 

The problem of (maximum) controlled invariant set computation for discrete-time systems 
has been a topic of active research for more than four decades. The central tool in this 
effort has been the contractive algorithm of [7] and its expansive counterpart [IH] . For an 
exhaustive survey and historical remarks see the survey [H] and the book |!3j . 

Both algorithms, although conceptually applicable to any nonlinear system, have been 
predominantly applied in a linear setting where they boil down to a sequence of linear 
programs and polyhedral projections. Finite termination of this sequence is a subtle prob- 
lem and sharp results are available only in the uncontrolled setting where no projections 
are required [17J; for discussion of finite-termination in the controlled case see [H] . The 
contractive and expansive algorithms were combined in [T5] to design an algorithm termi- 
nating in a finite number of iterations and outputting an e-accurate inner approximation 
of the MCI set (with the accuracy measured by the Hausdorff distance). Another line 
of research, culminating in |41j . exploits the linearity of the system dynamics in a more 
systematic way and approximates the maximum (or minimum) robust controlled invariant 
set by the Minkowski sum of a parametrized family of sets. Very recently, in continu- 
ous time, [31] developed a parallel algorithm for ellipsoidal approximations of the robust 
MCI set scalable to very high dimensions. Computation of low-complexity polyhedral 
controlled invariant sets was investigated in [TT] and [T2"] . 

In the nonlinear case, a common practice is to exploit the tight connection between in- 
variance and Lyapunov functions and seek invariant sets as sub-level sets of a (control) 
Lyapunov function; see, e.g., [T5l 150] and references therein for recent theoretical devel- 
opments on the related problem of region of attraction computation and, e.g., [33] for 
practical applications of these techniques. This, however, typically leads to non-convex 
bilinear optimization problems which are notoriously hard to solve. Therefore, one of- 
ten has to resort to ad-hoc analysis of the specific system at hand, which is typically 
tractable only in small dimensions; see (15J 06] for concrete examples. Related in spirit 
is the localization technique of |26j for discrete-time uncontrolled systems, also requiring 
considerable effort in analysing the system. 

Recently, a general approach using a hierarchy of finite-dimensional linear programs (LPs) 
was used in [6] to design a controller ensuring invariance of a given candidate polyhedral 
set. In our opinion, although being the current state of the art, this work still suffers 
from the following drawbacks: 1) the sets obtained are convex polytopes (not general 
semi-algebraic sets, a fact particularly limiting in the nonlinear case where nonconvex 
MCI sets are common); 2) the geometry of the candidate polytopic set must be given a 
priori; 3) there are no convergence guarantees to the MCI set. In this paper, we explicitly 
address all these points. 

Building upon our previous work [20] on the computation of the region of attraction (ROA) 
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for polynomial control systems, in this paper we characterize the maximum controlled 
invariant (MCI) set for discrete- as well as continuous-time polynomial systems as the 
solution to an infinite-dimensional LP problem in the cone of nonnegative measures. The 
dual of this problem is an infinite-dimensional LP in the space of continuous functions. 
Finite-dimensional relaxations of the primal LP and finite-dimensional approximations of 
the dual LP turn out to be semidefinite programs (SDPs) also related by duality. The 
primal relaxations lead to a truncated moment problem while the dual approximations to 
a sum-of-squares (SOS) problem. Super-level sets of one of the polynomials appearing in 
the dual SOS problem then provide outer approximations to the MCI set with guaranteed 
convergence as the degree of the polynomial tends to infinity. 

The main mathematical tool we use are the so-called occupation measures which allow 
us to study the time evolution of the whole ensemble of initial conditions (described by a 
measure) rather than studying trajectories associated to each initial condition separately. 
The use of measures to study dynamical systems has a very long tradition: see [13] for 
probably the first systematic treatment]^] for purely discrete-time treatment see (23J Chap- 
ter 6]. To the best of the authors' knowledge our paper is the first one to use occupation 
measures for MCI set (approximate) computation. The MCI set was previously charac- 
terized using occupation measures in [TO], but there the characterization is rather indirect 
and not straightforwardly amenable to computation. Apart from the authors' work |20j . 
the related problem of region of attraction computation was tackled using measures in |51j . 
There, however, a very different approach was taken, not using occupation measures but 
rather analyzing convergence via discretization of the state-space and propagating the 
initial distribution by means of a discretized transfer operator. Here, instead, we employ 
the (discounted) occupation measure which captures the behaviour of the trajectories 
emanating from the initial distribution over the infinite time horizon. As a result, our 
approach requires no discretization and, contrary to |5TJ, provides true guarantees (not 
in an "almost-everywhere" or "coarse" sense) and, more importantly, is applicable in a 
controlled setting. Closely related to the occupation measures used here is the Rantzer's 
density [52] which was used in [40] to assess the stability of attractor sets of uncontrolled 
nonlinear systems. The approach, however, does not immediately yield approximations 
of the MCI set (or the region of attraction) and applies to uncontrolled systems only. 

Similar in spirit to our approach, from the dual viewpoint of optimization over functions, 
are the Hamilton- Jacobi approaches (e.g., [36, [37]). However, contrary to these methods, 
our approach does not require state-space discretization and comes with convergence 
guarantees. 

The contribution of our paper with respect to previous work on the topic can be summa- 
rized as follows: 

• we deal with fully general continuous-time and discrete-time polynomial dynamics 
under semi-algebraic state and control constraints; 

• our approximated MCI set is described by (the intersections of) polynomial super- 
level sets, including more restrictive classes (e.g. polytopes, ellipsoids, etc.); 

1 ln [43 , J. E. Rubio used Young measures [JS] rather than occupation measures, but the basic idea 
of "linearizing" a nonlinear problem by going into an infinite-dimensional space of measures is the same. 
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• we provide a convex infinite-dimensional LP characterization of the MCI set; 

• we describe a hierarchy of convex finite-dimensional SDPs to solve the LP with 
convergence guarantees; 

• our approach is simple and readily applicable in the sense that the approximations 
are the result of a single SDP with no additional data required apart from the 
problem description. 

The contribution with respect to our previous work [SU] can be summarised as follows: 

• in [20 J we compute the ROA, which is a related although different object: it is the 
set of all of initial conditions that can be steered to a given target set while satisfying 
state and control constraints. In particular, the MCI set differs from the ROA in 
the sense that we do not try to hit any target set at a given time but rather try to 
keep the state within a given set forever. Therefore we had to adapt our technique 
to deal explicitly with invariance; 

• in [20J we dealt with continuous-time systems only, whereas we can cope, with 
minor modifications, with discrete-time systems as well; we choose to describe both 
the continuous-time and discrete-time setups in parallel precisely to underline these 
common features; 

• in [20J we considered only a finite time-horizon, whereas here we show how to cope, 
with the help of discounting, with an infinite horizon. This brought additional 
technical issues not encountered in finite time. 



What can be considered a drawback of our approach is the fact that the approximations 
to the MCI set we obtain are from the outside and therefore not invariant. However, 
accurate outer approximations provide important information as to the performance lim- 
itations of the control system and are of practical interest, e.g., in collision avoidance. 
Therefore we believe that our work bears both theoretical and practical value, and natu- 
rally complements existing inner-approximation techniques. 

The paper is organised as follows. The problem to be solved is described in Section [2] 
Occupation measures are introduced in Section [3] The infinite-dimensional primal and 
dual LPs are described in Sections|4]and[5j respectively. The finite-dimensional relaxations 
with convergence results are presented in Section [6} Numerical examples are in Section [7j 
A reader interested only in the semialgebraic outer approximations of the MCI set can 
consult directly the infinite-dimensional dual LPs ^ and ^ and their finite-dimensional 
approximations (11) and (13) in discrete and continuous time, respectively. 



1.1 Notation 

Measures are understood as signed Borel measures on a Euclidean space, i.e., as countably 
additive maps from the Borel sets to the real numbers. From now on all subsets of a 
Euclidean space we refer to are automatically understood as Borel. The vector space of 
all signed Borel measures with its support contained in a set X is denoted by M(X). 
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The support (i.e., the smallest closed set whose complement has a zero measure) of a 
measure « is denoted by spt u. The space of continuous functions on X is denoted by 
C(X) and likewise the space of once continuously differentiable functions is C l (X). The 
indicator function of a set X (i.e., a function equal to one on X and zero otherwise) is 
denoted by Ix{')- The symbol A denotes the n-dimensional Lebesgue measure (i.e., the 
standard n-dimensional volume). The integral of a function v with respect to a measure 
u over a set X is denoted by f x v(x) d/i(x). Sometimes for conciseness we use the shorter 
notation f vd/i omitting the integration variable and also the set over which we integrate 
if they are obvious from the context. The ring of polynomials in (possibly vector) variables 
Xi,. . . ,x n is denoted by WL[xi, . . . ,x n ]. 

2 Problem statement 

The approach is developed in parallel for discrete and continuous time. 

2.1 Discrete time 

Consider the discrete-time control system 

x t+ i = f(x t , u t ), x t eX, u t eU, t e {0, l, ...} (l) 

with a given polynomial vector field / with entries f\ G M.[x, u], i = 1, . . . ,n, and given 
compact basic semialgebraic state and input constraints 

x t EX :={xeR n : g Xi (x)>0,i = l,2,...,n x }, 
u t eU :={uer : gUi (u) > 0,i = 1,2,..., nu} 

with g Xi E R[x], 9Ui E R[u]. 

The maximum controlled invariant (MCI) set is defined as 

X I := jxo e X : 3 ({x t }Zi, MZi) s.t. x t +i = f{x t ,u t ), 

^eu, x t ex, vt g {0,1,...}}. 

A control sequence {u t }^ is called admissible if u t G U for all t 6 {0,1,...}. 

In words, the MCI set is the set of all initial states which can be kept inside the constraint 
set X ad infinitum using admissible control inputs. 

2.2 Continuous time 

Consider the relaxed continuous-time control system 

x(t) G conv f(x(t), U), x(t) G X, t G [0, oo), (2) 
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where conv denotes the convex hull, / is a polynomial vector field with entries G 
]R[x, u], i = 1, . . . ,n, and compact basic semialgebraic state and input constraint sets are 
defined by 

X :={xeR n : g Xi (x) > 0, i = 1, 2, . . . , n x }, 
U := {u G R m : g Ut {u) > 0, i = 1, 2, . . . , n v } 

with gxi G G M[u]. The meaning of the convex differential inclusion ^ is as 

follows: for all time t, the state velocity x(t) is constrained to the convex hull of the 
set f(x(t),U) := {f(x(t),u) : u G U} C M™. The connection of this convexified (or 
relaxed) control problem ^ and the classical control problem x = f(x, u) is the Filippov- 
Wazewski Theorem [5], which shows that the trajectories of x = f{x,u) are dense (in the 
supremum norm) in the set of trajectories of the convexified inclusiorj^] Therefore, 
from a practical point of view, there is little difference between the two formulations for 
the purposes of MCI set computation; see Section 3.2 and Appendices B and C of [20] for 
a detailed discussion on this subtle issue. The simplest assumption under which the MCI 
sets for both systems coincide is f(x, U) being convex for all x, which is in particular true 
for input-affine systems of the form x = f(x) + g(x)u with U convex. 

The maximum controlled invariant (MCI) set is defined as 

Xi := |x G X : 3 x(-) s.t. x(t) G conv/(x(t), U) a.e., x(t) G XWt G [0, oo) j, 

where x(-) is required to be absolutely continuous and a.e. stands for "almost everywhere" 
with respect to the Lebesgue measure on [0, oo). 

In words, the MCI set is the set of all initial states for which there exists a trajectory of 
the convexified inclusion ^ which remains in X ad infinitum. 

3 Occupation measures 

In this section we introduce the concept of occupation measures which is the centrepiece 
of our approach. 

3.1 Discrete time 

Given a discount factor a G (0,1), an initial condition xq and an admissible control 
sequence {u^}^ such that the associated state sequence {x t \ xo }tZo remains in X for all 
time, we define the discounted occupation measure /z(- | x Q ) G M(X x U) as 

oo 

x ) u t\x ) (3) 

t=o 

for all sets Ac X and B C U. 

2 Note that the set conv f(x(t), U) is closed for every t since / is continuous and U compact; therefore 
there is no need to take closure of the convex hull in order to apply the Filippov-Wazewski theorem. 
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In words, the discounted occupation measure measures the (discounted) number of visits 
of the state-control pair trajectory 1 x ), u{- | x )) to subsets of X x U . The discounting 
in the definition of the occupation measure ensures that fi(A x B \ xq) is always finite; in 
fact we have /i(X x U \ xq) = (1 — a) -1 . 

Now suppose that the initial condition is not a single point but an initial measure^ /a £ 
M(X) and an admissible control sequence is associated to each initial condition from the 
support of fio in such a way that the corresponding state sequence remains in X. Then 
we define the average discounted occupation measure /i £ M(X x U) as 

fi(AxB):= / n(A x B\x ) djj (x ). 
Jx 



The average discounted occupation measure measures the discounted average number of 
visits in subsets of X x U of trajectories starting from the initial distribution /io- 

Now we derive an equation linking the measures /io and /x. This equation will play a key 
role in subsequent development and in a sense replaces the dynamics equation Q. To 
derive this equation fix an initial condition xo <E X and a control sequence {u t \ xo }^ such 
that the associated state sequence {xt\x }t^o stays i n X. Then for any v £ C(X) we have 



„ oo oo 

/ v(x) dfi(x, u I x Q ) = ^ atv ( x t\x ) = v(x \ XQ ) + « ^ a t v(x t+1 

JxxU t=0 t=0 

oo 

= v(x 0{xo ) + a^2 atv (f( x t\ X0 ,u t i X0 )) 



t=0 



= v(x \ Xo ) + a / v(f(x,u))d[i(x,u\x ). 

J XxU 

Integrating w.r.t. /i we arrive at the sought equation 

/ v(x) d[i(x,u) = / v (x) dfio(x) + a / v{f{x 1 u))diJ,{x,u) Vf £ C(X). (4) 

JXxU Jx JXxU 



Note that this is an infinite-dimensional linear equation in variables (/io,/i). 

The following crucial Lemma establishes the connection between the support of any initial 
measure /io solving Q and the MCI set Xj. 

Lemma 1 For any pair of measures (/io, jj) satisfying equation with spt /io C X and 
spt /i C U x X we have spt /io C Xj. 



Proof: A detailed proof is in Appendix A. □ 

3 The initial measure /io can be thought of as the probability distribution of the initial state, although 
we do not require the mass of /iq to be normalized to one. 
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3.2 Continuous time 



Given an initial condition xq and a trajectory x(- \ Xq) of the inclusion (J5]) that remains 
in X for all £ > 0, there exists an admissible time- varying measure- valued relaxed control 
^t(-|%) G M(U), VtiUlxo) = 1, such that 

x(t) = / f(x(t), u) du t (u\x ) 
Ju 

almost everywhere with respect to the Lebesgue measure on [0, oo). This follows from 
the definition of the convex hull (in fact, for each t, v t {- \ x$) can be taken to be a convex 
combination of finitely many Dirac measures). 

Then, given a discount factor > 0, we define the discounted occupation measure //(• | xq) G 
M(X x U) as 

POO P 

H(A xB\x ) := / / e~ l3t I AxB (x(t\x ),u)dis t (u\x )dt 
Jo Ju 

for all sets Ac X and B cU. 

In words, the discounted occupation measure measures the (discounted) time spent by 
the state-control pair trajectory (x(- \x ),is(- \ x )) in subsets of X X U. The discounting 
in the definition of the occupation measure ensures that fi(A x B | xq) is always finite; in 
fact we have fi(X x U \ xq) = f3^ 1 . 

Now suppose that the initial condition is not a single point but an initial measure^ j2 e 
M(X) and a state trajectory that remains in X along with an admissible relaxed control 
is associated to each initial condition from the support of /j,q. Then we define the average 
discounted occupation measure fi e M(X x U) as 



fi(A x B) := / fx(A x B \ x ) dfx (x ). 
Jx 



Now we derive an equation linking the measures /j,q and \i. This equation will play a key 
role in subsequent development and in a sense replaces the dynamics equation (|2j). To 
derive the equation, fix an initial condition xq G X, a trajectory x(- \ Xq) that remains 
in X with an associated admissible relaxed control u t (- \ xq)- Then for any v G C 1 (X) 
integration by parts yields 



rad v ■ f(x, u) dfi(x, u \ xq) = j j e ^gr&dv ■ f(x(t \ xq), u) du t (u \ xq) dt 

d 



XxU JO JU 

oo 







e~ pt -r v{x{t\x ))dt 

iXv 



{3 e ^v(x(t \x )) dt — v(x(0\x )) 





P / v(x) dfx(x, u | x ) — v(x(0 | x )), 

J XxU 



4 The initial measure /io can be thought of as the probability distribution of the initial state, although 
we do not require the mass of /kq to be normalized to one. 
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where the boundary term at infinity vanishes due to discounting and the fact that X is 
bounded. Integrating with respect to /io then gives the sought equation 

/3 / v(x) dfx(x, u) = I v(x) dfio(x) + / gradf • f(x, u) dfi(x, u) Vv € C X (X). (5) 
JXxU Jx Jxxu 

Note that this is an infinite-dimensional linear equation in variables (/io,/t). 

The following crucial Lemma establishes the connection between the support of any initial 
measure satisfying ^ and the MCI set Xj. 

Lemma 2 For any pair of measures (/io, /•*) satisfying equation with spt /io C X and 
spt fi C U x X we have A(spt /i ) < A(Xj). 

Proof: A detailed proof is in Appendix B. □ 



4 Primal LP 

In this section we show how the MCI set computation problem can be cast as an infinite- 
dimensional LP problem in the cone of nonnegative measures. As in [20], the basic idea 
is to maximize the mass of the initial measure /io subject to the constraint that it be 
dominated by the Lebesgue measure, that is, /i < A. System dynamics is captured by 
the equations Q and ^ for discrete and continuous times, respectively; state and input 
constraints are expressed through constraints on the supports of the initial and occupation 
measure. The constraint that /U < A can be equivalently rewritten as /io + /io = A for 
some nonnegative slack measure /to G M(X). This constraint is in turn equivalent to 
J x w(x) dfi (x) + J x w(x) dfi (x) = j x w(x) d\(x) for all w G C(X). These considerations 
lead to the following primal LPs. 



4.1 Discrete time 

The primal LP in discrete time reads 

p* = sup /i (A) 

s.t. f v(x) dfi(x,u) = J v(x) dfio{x) + a J v (f(x, u)) dfi(x, u) WvEC(X) 
f w(x) dfj, (x) + f w(x) dfto(x) = f w(x) d\(x) Ww G C(X) 

^ > 0, /to > 0, /i > 
spt /i C X x U, spt fiQ C X, spt /i C X, 

(6) 

where the supremum is over the vector of measures (/i, /i ,/to) G M(X x U) x M(X) x 
M(X). 

This is an infinite-dimensional LP in the cone of nonnegative Borel measures. The fol- 
lowing Lemma, which is our main theoretical result, relates an optimal solution of this 
LP to the MCI set X/. 
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Theorem 1 The optimal value of LP problem is equal to the volume of the MCI set 
Xj, that is, p* = X(Xj). Moreover, the supremum is attained by the restriction of the 
Lebesgue measure to the MCI set Xj. 

Proof: The proof follows from Lemma[I]by the same arguments as Theorem 1 in J2DJ. By 
definition of the MCI set Xi, for any initial condition xq G Xj there exists an admissible 
control sequence such that the associated state sequence remains in X. Therefore for any 
initial measure /io < A with spt \xq C Xj there exist a discounted occupation measure \i 
with spt 11 C X x U and a slack measure /to with spt /to C X such that the constraints of 
problem ^ are satisfied. One such measure /xq is the restriction of the Lebesgue measure 
to Xj, and therefore p* > X(Xj). The fact p* < X(Xj) follows from Lemma [T] □ 



4.2 Continuous time 

The primal LP in continuous time reads 
p* = sup fi {X) 

s.t. f v(x) dfi(x,u) = J v(x) d[/,o(x) + f grad v ■ f(x, u) dfi(x, u) \/vEC l (X) 
f w(x) dfi (x) + f w(x) dfio(x) = f w(x) d\(x) Ww G C(X) 

fJ> > 0, (J,o > 0, /t > 
spt /i C X x U, spt /i c X, spt /t C X, 

(7) 

where the infimum is over the vector of measures (//, /j,q, /to) G M(X xU)x M(X) x M(X). 

This is an infinite-dimensional LP in the cone of nonnegative Borel measures. The fol- 
lowing Lemma, which is our main theoretical result, relates an optimal solution of this 
LP to the MCI set X/. 

Theorem 2 The optimal value of LP problem |?]) is equal to the volume of the MCI set 
Xj, that is, p* = X(Xj). Moreover, the supremum is attained by the restriction of the 
Lebesgue measure to the MCI set Xj. 

Proof: The fact that no equal to the restriction of the Lebesgue measure to Xj is feasible 
in (J7|) (and therefore p* > X(Xj)) follows by the same arguments as in discrete time. The 
fact that p* < X(Xj) follows from Lemma [2j □ 



5 Dual LP 

In this section we derive LPs dual to the primal LPs ^ and ([7]). Since the primal LPs 
are in the space of measures, the dual LPs will be on the space of continuous functions. 
Super-level sets of feasible solutions to these LPs then provide outer approximations to the 
MCI sets, both in discrete and in continuous time. Both duals can be derived by standard 
infinite-dimensional LP duality theory; see [20] for a derivation in a similar setting or [3] 
for a general theory of infinite-dimensional linear programming. 
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5.1 Discrete time 



The dual LP in discrete time reads 

d* = inf / w(x) d\(x) 
Jx 

s.t. av(f(x,u)) < v(x), V(x,u) eX xU (8) 
w(x) > v(x) + 1, V x G X 
w(x) > 0, Vx G X, 

where the infimum is over the pair of functions (v , w) G C(X) x C(X). 

The following key observation shows that the unit super-level set of any function w feasible 
in (|8f provides an outer-approximation to Xj. 

Lemma 3 Any feasible solution to problem |#[) satisfies v > and w > 1 on Xj. 

Proof: Given any x G X/ there exists a sequence {wt}^ , u t G U, such that Xf G X for 
all t. The first constraint of problem ^ is equivalent to av(x t+ i) < v(x t ), t G {0, 1, . . .}. 
By iterating this inequality we get 

v(x ) > a t v(xt) — > as t — > oo 

since x t G X and X is bounded. Therefore v(xo) > and w(xq) > 1 for all xq G Xj. □ 
The following theorem is instrumental in proving the convergence results of Section [6j 

Theorem 3 There is no duality gap between primal LP problems on measures and 
dual LP problem (^) on functions in the sense that p* = d* . 

Proof: Follows by the same arguments as Theorem 2 in [2D] using standard infinite- 
dimensional LP duality theory (see, e.g., [3]) and the fact that the feasible set of the 
primal LP is nonempty and bounded in the metric inducing the weak-* topology on 
M(X) x M(X x U) x M(X). To see non-emptiness, notice that the vector of measures 
(Ho, /-<■, Ao) = (0, 0, A) is trivially feasible. To see the boundedness, it suffices to evaluate the 
equality constraints of ^ for v(x) = w(x) = 1. This gives /io(X) + fto(X) = A(X) < oo 
and /i(X) = /i (X)/(l — a), which, since a G (0,1) and all measures are nonnegative, 
proves the assertion. □ 



5.2 Continuous time 

The dual LP in continuous time reads 

d* = inf J w(x) d\(x 
x 



s.t. grad-u • f(x, u) < (3v(x), V(x,w)gXx[/ (9) 
w(x) > v(x) + 1, Vx G X 

w(x)>0, VxgX, 
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where the infimum is over the pair of functions (v,w) G C l (X) x C(X). 

The following key observation shows that the unit super-level set of any function w feasible 
in ([9]) provides an outer-approximation to Xj. 

Lemma 4 Any feasible solution to problem |$|) satisfies v > and w > 1 on Xj. 

Proof: Given any x® G Xj there exists an admissible relaxed control function z^(-), 
Vt(U) = 1, such that x(t) G X for all t. For that x{t) we have J^t>(x(t)) = J^gradv • 
f(x(t),u)di , t (u) < fjj (3v (x(t)) dv t {u) = h , t {U)f3v(x{t)) = (3v(x(t)). Then by Gronwall's 
inequality v(x(t)) < e^v(xo), and consequently 

v(x ) > e~ /st v(x(t)) as t -> 00 

since x(t) G X and X is bounded. Therefore v(xq) > and w(xo) > 1 for all x$ G X/. □ 
The following theorem is instrumental in proving the convergence results of Section [6j 

Theorem 4 There is no duality gap between primal LP problems |?|) on measures and 
dual LP problem |Pp on functions in the sense that p* = d* . 

Proof: Follows by the same arguments as Theorem 2 in [20] using standard infinite- 
dimensional LP duality theory (see, e.g., and the fact that the feasible set of the 
primal LP is nonempty and bounded in the metric inducing the weak-* topology on 
M(X) x M(X x U) x M(X). To see non-emptiness, notice that the vector of measures 
(/i , yU 5 Ao) = (0, 0, A) is trivially feasible. To see the boundedness, it suffices to evaluate the 
equality constraints of ^ for v(x) = w(x) = 1. This gives ^o(X) + /to(X) = A(X) < 00 
and /i(X) = /i (X)//3, which, since /3 > and all measures are nonnegative, proves the 
assertion. □ 



6 LMI relaxations 

In this section we present finite dimensional relaxations of the infinite-dimensional LPs. 
Both in continuous and discrete time, the relaxations of the primal LPs lead to a truncated 
moment problem which translates to a semidefinite program (SDP) that can be solved 
by freely available software, e.g., SeDuMi [3H] or SDPA |4T] . Dual to the primal SDP 
relaxation is a sum-of-squares (SOS) problem that again translates to an SDP problem. 
The following discussion closely follows the one in |28j . 

We only highlight the main ideas behind the derivation of the finite-dimensional relax- 
ations. The reader is referred to [SHI Section 5] or to the comprehensive reference [32] 
for details. First, since the supports of all measures feasible in ^ and ([7]) are compact, 
these measures are uniquely determined by their moments, i.e., by integrals of all mono- 
mials (which is a sequence of real numbers when indexed in, e.g., the canonical monomial 
basis). Therefore, it suffices to restrict the test functions w(x) and v(x) in @ and @ to 
all monomials, reducing the linear equality constraints on measures /i , /i and fi of (|6| 
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and ([7]) to linear equality constraints on their moments. Next, by the Putinar Positivstel- 
lensatz (see [32J EH]), the constraint that the support of a measure is included in a given 
compact basic semialgebraic set is equivalent to the feasibility of an infinite sequence of 
LMIs involving the so-called moment and localizing matrices, which are linear in the co- 
efficients of the moment sequence. By truncating the moment sequence and taking only 
the moments corresponding to monomials of total degree less than or equal to 2k, where 
k G {1, 2, . . .} is the relaxation order, we obtain a necessary condition for this truncated 
moment sequence to be the first part of a moment sequence corresponding to a measure 
with the desired support. 

In what follows, M. k [-] denotes the vector space of real multivariate polynomials of total 
degree less than or equal to k. Furthermore, throughout the rest of this section we make 
the following standard standing assumption: 

Assumption 1 One of the polynomials modeling the sets X resp. U is equal to g Xi (x) — 
R\ — \\x\\1 resp. gui(u) = R\j — IMI2 R X , Ru sufficiently large constants. 

This assumption is completely without loss of generality since redundant ball constraints 
can be always added to the description of the compact sets X and U. 



6.1 Discrete time 



The primal relaxation of order k in discrete time reads 



p* k = max (|/o)o 

s.t. A k (y,y ,y ) 



M k (y)hO, M k „ dxi (g Xi ,y)hO, i = l,2,...,n x 

M k - du .{g Vi ,y) >z 0, i = 1,2, ...,n v 

M k (y o )h0, M Wx .( te ,2/o)hO, i = 1,2,. . . ,n x 

M k (y )yQ, M k _ dxt (g Xi ,y )yQ, % = 1, 2, . . . , n x , 



(10) 



where the notation >z stands for positive semidefinite and the minimum is over moment 
sequences (y, yo, yo) truncated to degree 2k corresponding to measures /i, /io and /to in (j6]). 
The linear equality constraint captures the two linear equality constraints of ^ with 
v (t, x) G K2fc[ij^] and w(x) G M 2 fc[a;] being monomials of total degree less than or equal to 
2k. The matrices M k (-) are the moment and localizing matrices, following the notations 



of [32] or [2D]. In problem (10), a linear objective is minimized subject to linear equality 



constraints and LMI constraints; therefore problem (10) is a semidefinite program (SDP). 
The dual relaxation of order k in discrete time reads 
d* k = inf w'l 

s.t. v(x) - av(f(x, u)) = q (x, u) + Y%=i Qi(x, u)g Xi (x) + J™=i n(x, u)gui{u) 
w(x) - v(x) - 1 = p (x) + Y!i=iPi(. x )9Xi(x) 
w(x) = s (x) + J27=i Si(x)g Xi (x), 

(11) 

where I is the vector of Lebesgue moments over X indexed in the same basis in which 
the polynomial w(x) with coefficients w is expressed. The minimum is over polynomials 
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~&2k[x] and w G M2fc[^], and polynomial sum-of-squares q^, Pi, s$, i — 1, . . . ,n x and 
1, . . . , njj, of appropriate degrees. In problem (11), a linear objective function is 
minimized subject to sum-of-squares (SOS) constraints; therefore problem (11) is an SOS 
problem which can be readily cast as an SDP (see, e.g., [32]). 



6.2 Continuous time 



The primal relaxation of order k in continuous time reads 



Pi 



max 
s.t. 



0o)o 

A k {y,y ,yo) 
M k (y) h 0, 



M k (y ) 
M h (y ) 



ho, 
h 0, 



M k _ dxi (g Xi ,y) hO, 
M k - dui (g Ui ,y) y 0, 
M k _ dxi (g Xi ,y ) h 0, 
M k _ dxi (g Xi ,y ) h 0, 



= 1,2,. 
= 1,2, 
= 1,2,. 
1,2,. 



.,n v 
.,n x 



:i2) 



where the notation >; stands for positive semidefmite and the minimum is over moment 
sequences (y, yo, yo) truncated to degree 2k corresponding to measures /i, /io and /to in (|7|. 
The linear equality constraint captures the two linear equality constraints of ^ with 
v(t, x) G IR-2fc [t, a?] and w(x) G M2fc[^] being monomials of total degree less than or equal to 
2k. The matrices M k (-) are the moment and localizing matrices, following the notations 



of [32] or [20]. In problem (12), a linear objective is minimized subject to linear equality 



constraints and LMI constraints; therefore problem (12) is a semidefinite program (SDP). 
The dual relaxation of order k in continuous time reads 



k 



inf w'l 
s.t. /3v(x) - 

w(x) — v(x) 

w(x) = Sq(x) 



radv- f(x,u) = q (x,u) + ^= 1 qi(x,u)g Xi (x) + ^=i r i( x > u )9u i (u 



[X 



E 



II X 

1=1 



{x)gxi(x), 



(13) 



where I is the vector of Lebesgue moments over X indexed in the same basis in which 
the polynomial w(x) with coefficients w is expressed. The minimum is over polynomials 
v(x) G M2fc[a^] and w G B^fc^], and polynomial sum-of-squares qi, pi, Si, % — 1, . . . , n x and 
rj, i = 1, . . . , njj, of appropriate degrees. In problem (13), a linear objective function is 



minimized subject to sum-of-squares (SOS) constraints; therefore problem (13) is an SOS 
problem which can be readily cast as an SDP (see, e.g., [3"2"]). 



6.3 Convergence results 



In this section we state several convergence results for the finite dimensional relaxations 



resp. approximations (10), (12) resp. (11), (13). Let w k and v k denote an optimal solution 
to the k th dual SDP approximation (11) or (13), and define 



X Ik :={xeX : v k (x) > 0}. 
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Then, in view of Lemmata [3] and |4j we know that w k over-approximates the indicator 
function of the MCI set Xj on X, i.e., w k > Ixi on X, and that the sets X Ik approximate 
from the outside the MCI set Xj, i.e., Xi k D X/. In the sequel we prove the following: 

• The optimal values of the finite-dimensional primal and dual problems p k and d* k 
coincide and converge to the optimal values of the infinite dimensional primal and 
dual LPs p* and d* which also coincide (in view of Theorems [3] and [1]) and are equal 
to the volume of the MCI set. 

• The sequence of functions w k converges on X from above to the indicator function 
of the MCI set in L\ norm. In addition, the running minimum minj<fc Wi converges 
on X from above to the indicator function of the MCI set set in L\ norm and almost 
uniformly. 

• The sequence of sets X Ik converges to the MCI set Xj in the sense that the volume 
discrepancy tends to zero, i.e., lim^oo \(X Ik \ Xj) = 0. 

The proofs of the results follow very similar reasoning as analogous results on region of 
attraction approximations in [201 Section 6]. 



Lemma 5 There is no duality gap between primal LMI problems (10 and 12) and dual 



L MI problems (11 and 13), i.e. p* k = d* k 



Proof: The argument closely follows the one in [20, Theorem 4] and therefore we only 
outline the key points of the proof. To prove the absence of duality gap, it is sufficient 



to show that the feasible sets of the primal SDPs (10) and (|12j) are non-empty and 
compact. The result then follows by standard SDP duality theory (see [20| Theorem 4] 
for a detailed argument). The non-emptiness follows trivially since the vector of measures 
(/i ,/i, /i) = (0, 0, A) is feasible in the primal infinite-dimensional LPs ^ and ^ and 
therefore the truncated moment sequences corresponding to these measures are feasible 



in the primal SDP relaxations (10) and (12). To see the compactness observe that the first 
components (i.e., masses) of the truncated moment vectors y®, y and y are bounded. This 
follows by evaluating the equality constraints of (|6| and ^ for w(x) = v(x) = 1. Indeed, 
in discrete-time we get (y) = (y ) Q /(l — a) and in continuous-time we get (y) = (yo)o/P) 
in addition, in both cases we have (?/o)o + (yo)o = A(X) < oo and therefore the first 
components are indeed bounded (since they are trivially bounded from below, in fact 
nonnegative, due to the constraints on moment matrices). Boundedness of the even 
components of each truncated moment vector then follows from the structure of the 
localizing matrices corresponding to the functions from Assumption [T] Boundedness of 
the entire truncated moment vectors then follows since the even moments appear on the 
diagonal of the positive semidefinite moment matrices. □ 

The following result shows the convergence of the optimal values of the relaxations to the 
optimal values of the infinite-dimensional LPs. 



Theorem 5 The sequence of infima of LMI problems (11) and (13) converges monoton- 
ically from above to the supremum of the LP problems^^ and (^J, i.e., d* < d* k+1 < d* k 



15 



and limfc^oo d* k 



d* . Similarly, the sequence of maxima of LMI problems |7o[ ) and (Hfy 
converges monotonically from above to the maximum of the LP problems ^ and (M), i.e., 
P* < P* k+1 < Pi and liiiifc^p* = p*. 



Proof: The monotonicity of the optimal values of the relaxations p* k resp. approximations 
d* k is evident form the structure of the feasible sets of the corresponding SDPs. The 
convergence of the primal relaxations pt to p* follows from the compactness of the feasible 



sets of the primal SDPs (10) and (12) (shown in the proof of Lemma pi) by standard 



arguments on the convergence of Lasserre's LMI hierarchy (see, e.g., (32j). The converge 
of the optimal value of the dual approximations d* k to d* then follows from Lemma [5j □ 

The next theorem shows functional convergence from above to the indicator function of 
the MCI set. 



Theorem 6 Let Wk G K^fcM denote the w-component of a solution to the dual LMI 
problems (11) or (IS) and let Wk{x) = minj<fc Wi(x). Then Wk converges from above to 
Ixj in L 1 norm and Wk converges from above to Ix, in L 1 norm and almost uniformly. 



Proof: The convergence in L\ norm follows immediately from Theorem [5] and from the 
fact that Wk > lx 1 by Lemmata [3] and |4j The convergence of the running minima follows 
from the fact that there exists a subsequence of {u>fc}^ which converges almost uniformly 
(by, e.g., [i Theorems 2.5.2 and 2.5.3]). □ 

Our last theorem shows a set-wise convergence of the outer-approximations to the MCI 
set. 

Theorem 7 Let (vk,Wk) € M 2 fc[x] x l^fcM denote an optimal solution to the dual LMI 



problem (11) or (13) and let X Ik := {x G M. n : v k (x) > 0}. Then Xi C X Ik , 



lim X(X Ik \X I ) = and X(n^ =1 X Ik \ X/) = 0. 

Proof: From Lemmata [3] or [4] we have Xj k D Xi and Wk > Ix^ therefore, since w > v + 1 
and w > on X, we have Wk > Ix Ik > Ix, and {x : Wk(x) > 1} D X Ik D X . From 
Theorem [6j we have Wk —> Ixj i n L 1 norm on X. Consequently, 



X(Xj) = / Ix, dX = lim / WkdX> lim / Ix Ik dX 
Jx k ^°° Jx k ^°° Jx 



lim X(Xjk) > lim A(n*L 1 X /i ) = A(n^ =1 X /fc ). 

k— >oo k — ^oo 



But since Xj C Xi k for all k, we must have 



lim X(X Ik ) = X(Xj) and A(n^ =1 X /fc ) = X(Xj) 

k—^oo 



and the theorem follows. 



□ 
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7 Numerical examples 



In this section we present numerical examples that illustrate our results. The primal 
SDP relaxations were modeled using Gloptipoly 3 (21] and the dual SOS problems using 
Yalmip P3J. The resulting SDP problems were solved using SeDuMi [38] (which, in the 
case of primal relaxations, also returns the dual solution providing the outer approxima- 
tions). For numerical computation (especially for higher relaxation orders), the problem 
data should be scaled such that the constraint sets are (within) unit boxes or unit balls; 
for ease of reproduction, most of the numerical problems shown are already scaled. On 
our problem class we observed only marginal sensitivity to the values of the discrete- and 
continuous-time discount factors a and and report results with a = 0.9 and (3 = 1 for 
all examples presented. 

For a discussion on the scalability of our approach and the performance of alternative 
SDP solvers see the Conclusion and the acrobot-on-a-cart example below. 

7.1 Discrete time 

7.1.1 Double integrator 

Consider the discrete-time double integrator: 

= X\ + 0.1X2 

x\ = X2 + 0.05u 

with the state constraint set X = [—1, l] 2 and input constraint set U = [—0.5,0.5]. The 
resulting of MCI set outer approximations of degree 8 and 12 are shown in Figure [T] the 
approximation is fairly tight for modest degrees. The true MCI set was computed using 
the standard algorithm based on polyhedral projections [9]. 

7.1.2 Cathala system 

Consider the Cathala system borrowed from [29J: 

X^ = X\ ~\~ X2 

x+ = -0.5952 + x 2 + x\. 

The chaotic attractor of this system is contained in the set X = [— 1.6, 1.6] 2 . MCI set 
outer approximations are shown in Figure [2] again, the approximations are relatively tight 
for small relaxation orders. The true MCI set was (approximately) computed by gridding. 

7.1.3 Julia sets 

Consider over z £ C, or equivalently over x £ IR 2 with z := X\ + ix 2 , the quadratic 
recurrence 

z + = z 2 + c 
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X 



X 



Figure 1: Discrete time double integrator - polynomial outer approximations (light gray) 
to the MCI set (dark gray) for degrees d G {8, 12}. 




Figure 2: Cathala system - polynomial outer approximations (light gray) to the MCI set (dark 
gray) for degrees d G {6, 10}. 

with c G C a given complex number and i the imaginary unit. The filled Julia set is the set 
of all initial conditions of the above recurrence for which the trajectories remain bounded. 
The shape of the Julia set depends strongly on the parameter c. If c lies inside the 
Mandelbrot set, then the Julia set is connected; otherwise the set is disconnected. In both 
cases the boundary of the set has a very complicated (in fact fractal) structure. Here we 
shall compute outer approximations of the filled Julia set intersected with the unit ball. To 
this end we set X = {x G 1R 2 : \\x\\ < 1}. Figure [3] shows outer approximations of degree 12 
for parameter values c = —0.7 + i0.2 (inside the Mandelbrot set) and c = —0.9 + z'0.2 
(outside the Mandelbrot set). The "true" filled Julia set was (approximately) obtained by 
randomly sampling initial conditions within the unit ball and iterating the recurrence for 
one hundred steps. Taking higher degree of the approximating polynomials does not give 
significant improvements due to our choice of the monomial basis to represent polynomials. 
An alternative basis (e.g. Chebyshev polynomials - see the related discussions in [22] and 
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c = 



-0.7 + i0.2 



c = 



-0.9 + z0.2 



Figure 3: Filled Julia set - polynomial outer approximation of degree 12 (light gray) and (an 
approximation of) the "true" set (dark grey) represented as an ensemble of initial conditions 
randomly sampled within the state-constraint set. The dashed line shows the boundary of the 
unit-ball state-constraint set. 

[20] ) would allow us to improve further the outer estimates and better capture the intricate 
structure of the filled Julia set's boundary. 

7.1.4 Henon map 

Consider the modified controlled Henon map 

x+ = 0.44 - 0.1x 3 - 4x 2 2 + 0-25u, 
x\ = Xi — 4xiX 2 , 
xt = x 2 , 

adapted from [33] with X = [— 1, l] 3 and U = [— M ma x, w max ]. We investigate two cases: 
uncontrolled (i.e., u max = 0) and controlled with M max = 1- Figure [5] shows outer approx- 
imations to the MCI set of degree eight for both settings and the "true" MCI set in the 
uncontrolled setting (approximately) obtained by random sampling of initial conditions 
inside the constraint set X. The outer approximations suggest that, as expected, allowing 
for control leads to a larger MCI set. 

7.2 Continuous time 
7.2.1 Double integrator 

Consider the continuous-time double integrator 

Xi = x 2 
x 2 = u, 
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Figure 4: Controlled Henon map - polynomial outer approximation of degree eight in the 
uncontrolled setting (darker red, smaller) and in the controlled setting (lighter red, larger). 
The (approximation of) the "true" set (black) in the uncontrolled setting is represented as an 
ensemble of initial conditions randomly sampled within the state-constraint set. 



with state constraint set X = [— 1, l] 2 and input constraint set U = [—1, 1]. The resulting 
MCI set outer approximations for degrees 8 and 12 are in Figure [5} The approximations 
are fairly tight even for relatively low relaxation orders. The true MCI set was (approxi- 
mately) computed as in Section 7.1.1 by methods of [9] after dense time discretization. 



7.2.2 Spider-web system 

As our second example we take the spider-web system from [1] given by equations 
x 1 = -0.15x[ + 200x^2 - 10.5x^x 2 - 807x^ + Ux\x\ + 600x 2 x^ - 3.5x^1 + 9x 



± 2 = -9x1 - 3.5x^X2 - 600x1x1 + \Ax\x\ + S07x\x A 2 - 10.5a; 2 x^ - 200x^1 - 0.15x 7 2 

with the constraint set X = [— 1, l] 2 . Here we exploit the fact that the system dynamics 
are captured by constraints on v only whereas w is merely over approximating v + 1, and 
the fact that outer approximations to the MCI set are given not only by {x : v{x) > 0} but 
also by {x : w(x) > 1}. Therefore, if low-complexity outer approximations are desired, it 



is reasonable to choose different degrees of v and w in (13) - high for v and lower for w - 
and use the set {x : w(x) > 1} as the outer approximation. That way, we expect to obtain 
relatively tight low-order approximations. This is confirmed by numerical results shown 
in Figure [6j The degree of v is equal to 16 for both figures, whereas deg w = 8 for the 
left figure and deg w = 16 for the right figure. We observe no significant loss in tightness 
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Figure 5: Continuous-time double integrator - polynomial outer approximations (light gray) to 
the MCI set (dark gray) for degrees d G {8, 14}. 

by choosing a smaller degree of w. The true MCI set was (approximately) computed by 
gridding. 




Figure 6: Spider-web system - polynomial outer approximations (light gray) to the MCI set 
(dark gray) for degrees degv = 16 and degw = 8 on the left and degw = 16 on the right. 



7.2.3 Acrobot on a cart 

As our last example we consider the acrobot on a cart system adapted from [23] , which is 
essentially a double pendulum on a cart where the inputs are the force acting on the cart 
and the torque in the middle joint of the double pendulum. The system is sketched in 
Figure [7| It is a sixth order system with with two control inputs; the dynamic equation 
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is given by 



x = 



x 5 

Xq 

M(x)~ l N(x,u) 



e 



6 



where 



a 2 cos x 2 



a 3 cos x 3 
a 5 cos(x 2 - x 3 ) 
a 6 



M(x) 



a 2 cos x 2 04 

03 cos x 3 a 5 cos(x 2 — X3) 



and 



N(x,u) 



U\ + a 2 X5 sin x 2 + a 3 x^ sin x 3 — (5 x 4 
•a 5 a;| sin(x 2 - x 3 ) + 5 2 x 6 + a 7 sinx 2 - x 5 (5x + 5 2 ) 
u 2 + a 5 sin(x 2 - x 3 )xl + 5 2 x 5 - 8 2 x 6 + a 8 sin x 3 



The states X\, x 2 , x 3 represent, respectively, the position of the cart (in meters), the 
angle of the lower rod and the angle of the upper rod of the double pendulum (both in 
radians); the states X4, £5 and x§ are then the corresponding velocities in meters per 
second for the cart and radians per second for the pendulum rods. The constants are 
given by a x = 0.85, a 2 = 0.2063, a 3 = 0.0688, a 4 = 0.0917, a 5 = 0.0344, a 6 = 0.0229, 
a 7 = 2.0233, a 8 = 0.6744, 5 = 0.3, 5\ = 0.1, 5 2 = 0.1. We are interested in computing 
the maximum controlled invariant subset of the state constraint set 



We investigate two cases. First, we consider the situation where only the middle joint is 
actuated and there is no force on the cart; therefore we impose the constraint (ui,u 2 ) G 
U = {0} x [—1, 1]. Second, we consider the situation where we can also exert a force on 
the cart; in this case we impose (ui,u 2 ) G U = [—1, 1] x [—1, 1]. Naturally, the MCI set 
for the second case is larger (or at least the same) as for the first case. This is confirmecj^] 
by outer approximations of degree four whose section for x\ = 0, 24 = 0, X5 = is shown 
in Figure [8] In order to compute the outer approximations we took a third order Taylor 
expansion of the non-polynomial dynamics even though exact treatment would be possible 
via a coordinate transformation leading to rational dynamics to which our methods can 
be readily extended; this extension is, however, not treated in this paper and therefore we 
opted for the simpler (and non-exact) approach using Taylor expansion. Before solving the 
problem we made a linear coordinate transform so that the state constraint set becomes 
the unit hypercube [— 1, l] 6 . 

This example, which is the largest of those considered in this paper, took 110 seconds 
to solvej^Jwith SeDuMi for d = 4; the corresponding time with the MOSEK SDP solver 
was 10 seconds. Using MOSEK we could also solve this example for d = 6 (in 420 
seconds) although there the solver converged to a solution with a rather poor accuracj|^] 
and therefore we do not report the results. 

5 There is no a priori guarantee on set- wise ordering of the outer approximations; what is guaranteed 
is the ordering of optimal values of the optimization problems (|12|) or (|13[) . 



6 A11 examples were run on an Apple iMac with 3.4 GHz Intel Corel?, 8 GB RAM and Mac OS X 
10.8.2. The time reported is the pure solver time, not including the Yalmip preprocessing time. 

7 Note that the MOSEK SDP solver is still being developed and its accuracy is likely to improve in 
the future. 



X = [-1,1] x [-7r/3,7r/3] x [-7r/3,7r/3] x [-0.5,0.5] x [-5,5] x [-5,5]. 
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Figure 7: Acrobot on a cart - sketch 




Figure 8: Acrobot on a cart - section of the polynomial outer approximations of degree four for 
(xi,X4, x$) = (0,0,0). Only the middle joint actuated - darker, smaller; middle joint and the 
cart actuated - lighter, larger. The states displayed x%, x% and xq are, respectively, the lower 
pendulum angle, the upper pendulum angle and the upper pendulum angular velocity. 

8 Conclusion 

We derived an infinite-dimensional convex characterization of the maximum controlled 
invariant (MCI) set, finite-dimensional approximations of (the dual of) which provide 
a converging sequence of semialgebraic outer-approximations to this set. The outer- 
approximations are the outcome of a single semidefinite program (SDP) with no addi- 
tional data required besides the problem description. Therefore the approach is readily 
applicable using freely available modeling tools such Gloptipoly 3 [2T] or YALMIP [31] 
with no hand-tuning involved. 

The cost to pay for this comfort is the relatively unfavourable scalability of the semidef- 
inite programs solved - the number of variables grows a.s 0((n + m) d ), where n and m 
are the state and control dimensions and d is the degree of the approximating polyno- 
mial. Therefore, in order for this approach to scale to medium dimensions (say, more 
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than m + n = 6) one either has to tradeoff accuracy by taking small d or go beyond 
the standard freely available solvers such as SeDuMi or SDPA. One possibility is par- 
allelization; for instance, the free parallel solver SDPARA [38] allows for the approach 
to scale to larger dimensions. Alternatively, one can utilize one of the (few) commercial 
SDP solvers; in particular, the recently released MOSEK SDP solver seems to show far 
superior performance on our problem class, and therefore this may allow the approach to 
scale to larger dimensions (see also the discussion following the acrobot-on-a-cart example 
in Section 7.2.3). Finally, one can resort to customized structure-exploiting solutions; this 
is a promising direction of future research currently investigated by the authors. At this 
point it should be emphasized that, to the best of the authors' knowledge, all of the ex- 
isting approaches providing approximations of similar quality experience similar or worse 
scalability properties. 

Other directions of future research include the extension of the presented approach to 
inner approximations of MCI sets, to stochastic systems and to uncertain systems. Par- 
tial results on the inner approximations for the related problem of region of attraction 
computation already exist [2H], albeit in uncontrolled setting only. 



Appendix A 

We start by embedding our problem in the setting of discrete-time Markov control pro- 
cesses; terminology and notation is borrowed from the classical reference [23J. Let us 
define a stochastic kernel on U given A as a map v{- | •) such that v[- \ x) is a prob- 
ability measure on U for all x G X and v[B | •) is a measurable function on X for all 
B C U. Any such stochastic kernel gives rise to a discrete-time Markov process when 
applied to system Q as a stationary randomized control policy (a policy which, given 
x, chooses the control action randomly based on the probability distribution u(- \ x), i.e., 
Prob(w G B | x) = v[B \ x) for all B C U). The transition kernel Q v {- | •) of this stationary 
Markov process is then given by 

Q v (A\x)= [ I A (f(x,u))dv(u\x) = Prob(x+ G A\x) VAcM n , 
Ju 

where x is the current state and x + the successor state. The t-step transition kernel is 
then defined by induction as 

Qt(A\x):= [ QiA^dQl^iylx), te{2,3,...} 

with Ql := Q u . Given an initial distribution /i , the distribution of the Markov chain at 
time t, fit, is given by 

(k(A)= / Ql(A\x) dfi Q (x) = Prob(x t G A). 

The joint distribution of state and control is then 

fj, t (A x B) — I u(B\x) djj t (x). 
J A 
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The discounted occupation measure associated to the Markov process is defined by 

oo 

H(A x B) =^o*/it(i4 x B). 



t=o 



Note that this relation reduces to ([3j) when ii t = o~(x t ,u t ) 

In order to prove Lemma [I] we need the following result that can be found in [23] . 

Lemma 6 For any pair of measures (/x , fJ>) satisfying equation O there exists a station- 
ary randomized control policy v{- \ x) such the Markov chain obtained by applying this 
control policy to the difference equation starting from initial distribution fi has the 
discounted occupation measure equal to fi. 

Proof: Disintegrate fi as dfi(x,u) = du(u \ x)dfi(x), where fx denotes the x-marginal 
of 11 and v is a stochastic kernel on U given X. According to the discussion preceding 
Lemma [6j applying v to (jl]) gives rise to a stationary discrete-time Markov process with 
the transition kernel Q u starting form the initial distribution /x . 

With this notation, equation Q can be equivalently rewritten as 



v (x) dfi(x) = / v (x) dfj, (x) + a / / v{y) dQ u (y\ x) dfx{x) (14) 

X JX Jx JR" 

for any measurable v(x) (derivation of equation Q did not depend on the continuity of v). 
Taking v{x) := Ia{%) we obtain 

Pl{A) = Li Q {A) + a Q v {A\x)dfi{x) VAcI. (15) 
Jx 



hand side of (15) we get 



Using relation (14) with v(x) := Q U (A \ x) to evaluate the integral w.r.t. fi on the right 



li(A) = fMo(A) + a / Q L ,(A\x) dfio(x) + a / Q u (A\x) djl(x). 
Jx Jx 

By iterating this procedure we obtain 

ft(A)=ti (A)+J2a i f Qi(A\x)dno(x) + a t+1 [ Ql +1 (A \ x) dfi(x), (16) 

m{A) o 

and taking the limit as t — > oo gives 

oo 

fi{A) = ^a%{A), 
t=o 

where the third term in (16) converges to zero because a G (0, 1), Q t l ^ rl {A \ x) < 1 and 
jl is a finite measure. Hence the x-marginal of the discounted occupation measure of the 
Markov chain coincides with the x- marginal of \x. 
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Finally, to establish equality of the whole measures observe that 

OO OO p p 

^a'/i ( (A x B) = ^a' / v{B\x) djx t {x) — / u(B \x) djx{x) — fx(A x B) 



□ 



Proof of Lemma [TJ Disintegrate lx to dji{x,u) = dv{u \ x)djx(x) as in the proof of 
Lemma [6} Then for any x G S : = spt jx we have 

I s (f(x,u)) v(u\x) = 1. 



This relation says that the support of jx is invariant under v and follows from Lemma [6j 
from the definition of the occupation measure fx, from the definition of the support and 
from the fact that z/(- \x) is a probability measure for all x. 

Define an admissible stationary deterministic control policy by taking any measurable 
selection u(x) G spt v(- \ x) C U. Define further the sequence of probability measures 

, . . u(B 1/n (u(x)) n A\x) 
Vn{A\x)= 1/n Vne 1,2,... , AcU, 

is(B 1/n (u{x)) nu\x) 

where By n (u(x)) is a closed ball of radius 1/n centered at u(x). Then v n {-\x) converges 
weakly-* (or weakly or narrowly) to S u m and 



I s {f{x,u))v n {u\x) = 1 Vn G {1,2,...}. 



v 



Therefore, 



l = limsup/ I s (f(x,u))u n (u\x) < / I s (f(x,u))5 u{x )(u) = I s (f(x,u(x))), 

n-^oa Jjj Ju 

where the inequality follows by the Portmanteau lemma since the set {u \ f(x,u) G 
S PI By n (u(x))} is closed for all x by continuity of /. Therefore in fact Is(f(x, u(x))) = 1 
and so f(x, u(x)) G spt jx for all x G spt jx. Therefore spt jx C X is invariant for the closed 
loop system xt+i = f(x t ,u(x t )), where u(x) is an admissible deterministic control policy. 
Therefore necessarily spt jx C Xj. Finally, from equation Q clearly spt llq C spt jx and so 
spt fXo C X/. □ 
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Lemma 7 For any pair of measures (/xq, lx) solving there exists a family of trajec- 
tories of the convexified inclusion (pjp starting from /xq such that the x-marginal of its 
discounted occupation measure is equal to the x-marginal of \x. 

Proof: The proof is based on fundamental results of j2] and [8] and on the compactifica- 
tion procedure discussed in [30] ■ 
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We begin by embedding the problem in a stochastic setting. To this end, define the 
extended state space E as the one-point compactification of R n , i.e., E = IR™U {A}, where 
A is the point compactifying R n . Define also the linear operator A : T>(A) — )■ C(E x Z7) 
by 

w I—)- Aw := gradw ■ /, 
where the domain of A, T>(A), is defined as 

V(A) :={w:E^R\we C 1 ^™), w(A) = 0, lim w(x) = 0, 

x— >A 

lim gradw • f(x,u) = Vm G (/}. 

In words, T>(A) is the space all continuously differentiable functions vanishing at infinity 
such that gradu> • / also vanishes at infinity for all u G U. Now consider the relaxed 
martingale problem [8]: find a stochastic process Y : [0, oo] x Q — > E defined on some 
filtered probability space (Q, J 7 , (J r t)t>o, P) an d a stochastic kernel u(- \ •) (stationary 
relaxed Markov control) on U given E such that 



• P(Y(0) e A) = hq(A) VAc£ 

• for all w e T>(A) the stochastic process 

w(Y(t))- [ [ Aw(Y(r),u)u(du\Y(T))dr (17) 
Jo Ju 

is an martingale (see, e.g., [23] for a definition). 

Observe that there exists a countable subset of T>(A) (e.g., all polynomials with rational 
coefficients attenuated near infinity) dense in T>(A) in the supremum norm. Next, T>(A) 
is clearly an algebra that separates points of E and Al = 0. Finally, since f(x,u) is 
polynomial and hence locally Lipschitz, the ODE x = f(x,u) has a solution on [0, oo) 
for any xq € E and any fixed u £ U in the sense that if there is a finite escape time 
t e , then we define x(t) = A for all t > t e . Each such solution satisfies the martingale 
relation (17) (with a trivial probability space). Therefore, A satisfies Conditions 1-3 of [8] 
and it follows from Theorem 2.2 and Corollary 2.2 therein that for any pair of measures 
satisfying the discounted Liouville's equation (|5]), there exists a solution to the above 
martingale problem whose discounted occupation measure is equal to //, that is, 

fi(A x5) = E{| e- pt I AxB (Y(t), u) u{du \ Y(t)) dt}, P(Y(0) e A) = ^(A), 

where E denotes the expectation w.r.t. the probability measure P. From the martingale 



property of (17) and the definition of A we get 

E{w(Y(t))}-E^J j gradw • f(Y(r),u) v{du\Y(r)) drj = E{y(0)}. 

Now let Ht denote the marginal distribution of Y(t) at time t; that is, 
m(A) := P(Y(t) eA)= E{I A (Y(t))} V A C X. 
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Then the above relation becomes 



/ w{x) dji t {x) — III grad w(x) ■ f(x, u) v{du \ x) dfx T (x) dr = / w(x) dfj,( ] (x), 
Jx Jo Jx Ju J 

where we have used Fubini's thorem to interchange the expectation operator and integra- 
tion w.r.t. time. Defining the relaxed vector field 



f(x)= / f(x, u) v(du | x) G conv f(x, U) 
Ju 

and rearranging we obtain 

w(x) dfx t (x) — / w(x) dpL (x) + / / gradw(:c) • f(x) d/j, T (x) dr, (18) 
x J Jo Jx 

where the equation holds for all w G C 1 (X) almost everywhere with respect to the 
Lebesgue measure on [0, oo). The Lemma then follows from Ambrosio's superposition 
principle |2, Theorem 3.2] using the same arguments as in the proof of Lemma 4 in [20] . 

□ 

Proof of Lemma [2j Suppose that a pair of measures (jUo>A*) satisfies ^ and that 
A(spt fiQ \Xj) > 0. From Lemma [7] there is a family of trajectories of (|2| starting from /j, 
with discounted occupation measure whose x-marginal coincides with the x-marginal of 
fi. However, this is a contradiction since no trajectory starting from spt fiQ \ Xj remains 
in X for all times and spt fi C X. Thus, A (spt /xq \ Xj) = and so A (spt /io) < X(Xi). □ 
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